---
title: "gender_typicality_analysis RR3"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE,message=FALSE)
setwd("/Users/martinnaunov/Desktop/Desktop - Martin’s MacBook Pro (2)/Sexuality and Gender Typicality 2024")
```


```{r, include=FALSE}
rm(list=ls())
```

```{r}
setwd("/Users/martinnaunov/Desktop/Desktop - Martin’s MacBook Pro (2)/Sexuality and Gender Typicality 2024")

gt_national <- read.csv("gt_national_working.csv")

gt_pssp <- read.csv("gt_pssp_working.csv")

```

```{r}
library(lmerTest)
```


```{r releveling reference categories}
gt_national$gt <- as.factor(gt_national$gt)

gt_national <- within(gt_national, gt <- relevel(gt, ref = "mass"))
gt_national$so<-as.factor(gt_national$so)
gt_national <- within(gt_national, so <- relevel(so, ref = "straight"))


gt_pssp$gt <- as.factor(gt_pssp$gt)

gt_pssp <- within(gt_pssp, gt <- relevel(gt, ref = "mass"))
gt_pssp$so<-as.factor(gt_pssp$so)
gt_pssp <- within(gt_pssp, so <- relevel(so, ref = "straight"))
```

```{r}
library(dplyr)
#Support National (Recoding incl Binary) 
gt_national$support_2<-recode(gt_national$support,"Very Unlikely" = 1, "Unlikely" = 2, "Somewhat Unlikely"=3, "Somewhat Likely"=4, "Likely"=5,"Very Likely"=6)

gt_national$Support <- as.factor(gt_national$support_2)

range01 <- function(x, ...){(x - min(x, ...)) / (max(x, ...) - min(x, ...))}


gt_national$support_binary<-dplyr::recode(gt_national$support,"Very Unlikely" = 0, "Unlikely" = 0, "Somewhat Unlikely"=0, "Somewhat Likely"=1, "Likely"=1,"Very Likely"=1)

#Support Student (Recoding incl Binary) 


gt_pssp$support_2 <- dplyr::recode(gt_pssp$support,
                           "Very Unlikely" = 1, "Unlikely" = 2, "Somewhat Unlikely"=3,
                           "Somewhat Likely"=4,"Likely"=5,"Very Likely"=6)

gt_pssp$Support <- as.factor(gt_pssp$support_2)


gt_pssp$support_binary<-recode(gt_pssp$support,"Very Unlikely" = 0, "Unlikely" = 0, "Somewhat Unlikely"=0, "Somewhat Likely"=1, "Likely"=1,"Very Likely"=1)

table(gt_pssp$support_binary)

```


# Feelings based on SO & GT

```{r Feelings based on SO & GT National}

library(lme4)

model_feeling_national_1 <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), 
                                           data = gt_national)
summary(model_feeling_national_1)



model_feeling_national_2 <- lmer(feeling ~ so*gt + (1 | respondent)+ (1 | candidate), 
                                           data = gt_national)
summary(model_feeling_national_2)


```


```{r Feelings based on SO & GT Student}
model_feeling_pssp_1 <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), 
                                           data = gt_pssp)
summary(model_feeling_pssp_1)


model_feeling_pssp_2 <- lmer(feeling ~ so*gt + (1 | respondent)+ (1 | candidate), 
                                           data = gt_pssp)
summary(model_feeling_pssp_2)

```


```{r}
#library(pbkrtest)

#KRmodcomp(model_feeling_national_1, model_feeling_national_2)

#Ftest    1.7452    2.0000 6553.8752         1  0.1747
```




# Support Based on SO & GT

```{r Support based on SO & GT National}
model_binary_national_1 <- glmer(support_binary ~ so + gt + (1|respondent)+(1|candidate),
                              data = gt_national, family=binomial(link='logit'))

summary(model_binary_national_1)



model_binary_national_2 <- glmer(support_binary ~ so*gt + (1|respondent)+(1|candidate),
                              data = gt_national, family=binomial(link='logit'))

summary(model_binary_national_2)


```



```{r Support based on SO & GT Student}

model_binary_pssp_1 <- glmer(support_binary ~ so + gt + (1|respondent)+(1|candidate),
                              data = gt_pssp, family=binomial(link='logit'))

summary(model_binary_pssp_1)


model_binary_pssp_2 <- glmer(support_binary ~ so*gt + (1|respondent)+(1|candidate),
                              data = gt_pssp, family=binomial(link='logit'))

summary(model_binary_pssp_2)

```


# Custom function to extract coefficients, confidence intervals, and p-values for lmer & glmer models
```{r}
# Custom function to extract coefficients, confidence intervals, and p-values for lmer models
extract_lmer_info <- function(model) {
  
  # Extract coefficients (fixed effects)
  coef_values <- summary(model)$coefficients[, "Estimate"]
  
  # Extract standard errors
  std_errors <- summary(model)$coefficients[, "Std. Error"]
  
  # Extract p-values (Satterthwaite's method in lmer models)
  p_values <- summary(model)$coefficients[, "Pr(>|t|)"]
  
  # Calculate 95% confidence intervals using confint()
  conf_intervals <- confint(model, parm = "beta_", level = 0.95)
  lower_ci <- conf_intervals[, 1]
  upper_ci <- conf_intervals[, 2]
  
  # Create a data frame with all the extracted information
  results <- data.frame(
    Term = rownames(summary(model)$coefficients),
    Estimate = coef_values,
    Lower_95CI = lower_ci,
    Upper_95CI = upper_ci,
    P_value = p_values
  )
  
  return(results)
}

# Custom function to extract coefficients, confidence intervals, and p-values for glmer models
extract_glmer_info <- function(model) {
  
  # Extract coefficients (log-odds)
  log_odds <- summary(model)$coefficients[, "Estimate"]
  
  # Extract standard errors
  std_errors <- summary(model)$coefficients[, "Std. Error"]
  
  # Extract p-values
  p_values <- summary(model)$coefficients[, "Pr(>|z|)"]
  
  # Calculate 95% confidence intervals for log-odds
  z_value <- 1.96
  lower_ci_log_odds <- log_odds - z_value * std_errors
  upper_ci_log_odds <- log_odds + z_value * std_errors
  
  # Exponentiate the log-odds to get odds ratios and confidence intervals
  odds_ratios <- exp(log_odds)
  lower_ci_odds_ratios <- exp(lower_ci_log_odds)
  upper_ci_odds_ratios <- exp(upper_ci_log_odds)
  
  # Create a data frame with all the extracted information
  results <- data.frame(
    Term = rownames(summary(model)$coefficients),
    Estimate = odds_ratios,  # Now using "Estimate" instead of "Odds_Ratio"
    Lower_95CI = lower_ci_odds_ratios,
    Upper_95CI = upper_ci_odds_ratios,
    P_value = p_values
  )
  
  return(results)
}


```

## Custom function to extract variances and group sizes from lmer/glmer models
```{r}
# Custom function to extract variances and group sizes from lmer/glmer models
extract_variances_lmer <- function(model) {
  # Extract variance components
  variance_components <- as.data.frame(VarCorr(model))
  
  # Extract the number of groups (respondents and candidates)
  num_respondents <- summary(model)$ngrps["respondent"]
  num_candidates <- summary(model)$ngrps["candidate"]
  
  # Extract respondent and candidate variances
  respondent_variance <- variance_components[variance_components$grp == "respondent", "vcov"]
  candidate_variance <- variance_components[variance_components$grp == "candidate", "vcov"]
  residual_variance <- variance_components[variance_components$grp == "Residual", "vcov"]
  
  # Create a dataframe with the extracted information
  variance_info <- data.frame(
    Model = deparse(substitute(model)),
    Respondent_Variance = respondent_variance,
    Candidate_Variance = candidate_variance,
    Residual_Variance = residual_variance,
    Num_Respondents = num_respondents,
    Num_Candidates = num_candidates
  )
  
  return(variance_info)
}

extract_variances_glmer <- function(model) {
  # Extract variance components
  variance_components <- as.data.frame(VarCorr(model))
  
  # Extract the number of groups (respondents and candidates)
  num_respondents <- summary(model)$ngrps["respondent"]
  num_candidates <- summary(model)$ngrps["candidate"]
  
  # Extract respondent and candidate variances
  respondent_variance <- variance_components[variance_components$grp == "respondent", "vcov"]
  candidate_variance <- variance_components[variance_components$grp == "candidate", "vcov"]
  
  # Create a dataframe with the extracted information
  variance_info <- data.frame(
    Model = deparse(substitute(model)),
    Respondent_Variance = respondent_variance,
    Candidate_Variance = candidate_variance,
    Num_Respondents = num_respondents,
    Num_Candidates = num_candidates
  )
  
  return(variance_info)
}

```


#Table 1: The Effect of Candidate Sexuality and Gender Presentation on Electoral Viability
```{r}
format_model <- function(estimate, lower_ci, upper_ci, p_value) {
  # Format the estimate (OR for binary models, Estimate for lmer models)
  estimate_str <- sprintf("%.2f", estimate)
  
  # Format the confidence intervals
  ci_str <- sprintf("(%.2f, %.2f)", lower_ci, upper_ci)
  
  # Add stars for p-value significance
  stars <- ifelse(p_value < 0.001, "***",
                  ifelse(p_value < 0.01, "**",
                         ifelse(p_value < 0.05, "*", "")))
  
  # Concatenate the estimate, stars, and confidence intervals
  formatted_result <- paste0(estimate_str, stars, "\n", ci_str)
  
  return(formatted_result)
}

```

```{r}
# Apply the functions to the glmer support models
results_national_1 <- extract_glmer_info(model_binary_national_1)
results_national_2 <- extract_glmer_info(model_binary_national_2)
results_pssp_1 <- extract_glmer_info(model_binary_pssp_1)
results_pssp_2 <- extract_glmer_info(model_binary_pssp_2)

# Apply the function to the lmer feeling models
results_feeling_national_1 <- extract_lmer_info(model_feeling_national_1)
results_feeling_national_2 <- extract_lmer_info(model_feeling_national_2)
results_feeling_pssp_1 <- extract_lmer_info(model_feeling_pssp_1)
results_feeling_pssp_2 <- extract_lmer_info(model_feeling_pssp_2)


# Add model names to each dataframe
results_national_1$Model <- "National Model 1 (OR)"
results_national_2$Model <- "National Model 2 (OR)"
results_pssp_1$Model <- "PSSP Model 1 (OR)"
results_pssp_2$Model <- "PSSP Model 2 (OR)"
results_feeling_national_1$Model <- "Feeling National Model 1"
results_feeling_national_2$Model <- "Feeling National Model 2"
results_feeling_pssp_1$Model <- "Feeling PSSP Model 1"
results_feeling_pssp_2$Model <- "Feeling PSSP Model 2"


# Combine all results into one dataframe
combined_results <- rbind(results_national_1, results_national_2, results_pssp_1, results_pssp_2,
                          results_feeling_national_1, results_feeling_national_2, 
                          results_feeling_pssp_1, results_feeling_pssp_2)

# Print combined results
print(combined_results)

```


```{r}

combined_results <- combined_results %>%
  mutate(formatted_result = format_model(Estimate, Lower_95CI, Upper_95CI, P_value))

# Print the modified dataframe
print(combined_results)

# Reshape the dataframe to a wide format with the models as columns
wide_results <- combined_results %>%
  select(Term, Model, formatted_result) %>%
  tidyr::pivot_wider(names_from = Model, values_from = formatted_result)


# Reorder the columns in the desired order
wide_results <- wide_results %>%
  select(Term, 
         "PSSP Model 1 (OR)", "PSSP Model 2 (OR)", 
         "Feeling PSSP Model 1", "Feeling PSSP Model 2", 
         "National Model 1 (OR)", "National Model 2 (OR)", 
         "Feeling National Model 1", "Feeling National Model 2")

# Print the reordered dataframe
print(wide_results)


```


```{r}
library(kableExtra)

# Create LaTeX table with the reordered columns
latex_table <- kable(wide_results, format = "latex", booktabs = TRUE, escape = FALSE,
                     col.names = c("Term", 
                                   "PSSP Model 1 (OR)", "PSSP Model 2 (OR)", 
                                   "Feeling PSSP Model 1", "Feeling PSSP Model 2", 
                                   "National Model 1 (OR)", "National Model 2 (OR)", 
                                   "Feeling National Model 1", "Feeling National Model 2")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down")) %>%
  add_header_above(c(" " = 1, "PSSP Models" = 4, "National Models" = 4))  # Custom header for the models

# Save as a LaTeX file
cat(latex_table, file = "table_1_so_gt_results.tex")

```

#Second Part of Table 1
```{r}
variance_pssp_1 <- extract_variances_glmer(model_binary_pssp_1)
variance_pssp_2 <- extract_variances_glmer(model_binary_pssp_2)
variance_feeling_pssp_1 <- extract_variances_lmer(model_feeling_pssp_1)
variance_feeling_pssp_2 <- extract_variances_lmer(model_feeling_pssp_2)
variance_national_1 <- extract_variances_glmer(model_binary_national_1)
variance_national_2 <- extract_variances_glmer(model_binary_national_2)
variance_feeling_national_1 <- extract_variances_lmer(model_feeling_national_1)
variance_feeling_national_2 <- extract_variances_lmer(model_feeling_national_2)


variance_pssp_1$Residual_Variance <- NA
variance_pssp_2$Residual_Variance <- NA
variance_national_1$Residual_Variance <- NA
variance_national_2$Residual_Variance <- NA

# Combine the results into a single dataframe
combined_variances <- rbind(variance_pssp_1, variance_pssp_2, variance_feeling_pssp_1, variance_feeling_pssp_2, variance_national_1, variance_national_2, variance_feeling_national_1, variance_feeling_national_2
)

# Print combined variances dataframe
print(combined_variances)

```


```{r}
library(tidyr)
library(dplyr)

# Reshape the combined_variances dataframe so that the models are columns and variances are rows
reshaped_variances <- combined_variances %>%
  pivot_longer(cols = c(Respondent_Variance, Candidate_Variance, Residual_Variance, Num_Respondents, Num_Candidates), 
               names_to = "Metric", 
               values_to = "Value") %>%
  pivot_wider(names_from = Model, values_from = Value)

# Print reshaped variances
print(reshaped_variances)

```

```{r}

reshaped_variances <- reshaped_variances %>%
  mutate(across(where(is.numeric), ~ sprintf("%.2f", .)))


# Generate the LaTeX table
latex_table_variances <- reshaped_variances %>%
  kable(format = "latex", booktabs = TRUE, escape = FALSE, col.names = c("Metric", "PSSP Model 1 (OR)", "PSSP Model 2 (OR)", 
                                                                        "Feeling PSSP Model 1", "Feeling PSSP Model 2", 
                                                                        "National Model 1 (OR)", "National Model 2 (OR)", 
                                                                        "Feeling National Model 1", "Feeling National Model 2")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down"))

# Save the LaTeX table to a file (if you want)
cat(latex_table_variances, file = "reshaped_variances_table.tex")


```


#Predicted Probability (Figure 4)

```{r}
# Load necessary libraries
library(ggplot2)
library(dplyr)
library(ggpattern)
library(cowplot)
library(lme4)
library(ggsignif)
library(ggtext)

# Fit the model using glmer
model_binary_national_1 <- glmer(support_binary ~ so + gt + (1|respondent) + (1 | candidate),
                                 data = gt_national, family = binomial(link = 'logit'))

# Calculate predicted probabilities and add to the dataset
gt_national$predicted_prob <- predict(model_binary_national_1, type = "response")

# Reorder levels of 'gt' (Gender Typicality) to be 'mass', 'neutral', 'fem'
gt_national$gt <- factor(gt_national$gt, levels = c("mass", "neutral", "fem"))


# Function to compute summary statistics with confidence intervals
compute_summary <- function(data, group_vars) {
  data %>%
    group_by(across(all_of(group_vars))) %>%
    summarize(
      mean_predicted_prob = mean(predicted_prob),
      se = sd(predicted_prob) / sqrt(n()),  # Standard error
      lower = mean_predicted_prob - 1.96 * se,  # 95% CI lower bound
      upper = mean_predicted_prob + 1.96 * se,  # 95% CI upper bound
      .groups = "drop"
    )
}


y_limits <- c(0, 0.80)

plot1_data <- compute_summary(gt_national, "so")

levels(plot1_data$so) <- c("Straight", "Gay")

plot1 <- ggplot(plot1_data, aes(x = so, y = mean_predicted_prob, pattern = so)) +
  geom_bar_pattern(
    stat = "identity", 
    color = "black", 
    width = 0.5,
    pattern_fill = "black",
    pattern_density = 0.1,
    pattern_spacing = 0.05,
    fill = "white"
    ) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, color = "black") +
  geom_signif(
    comparisons = list(c("Straight", "Gay")),
    annotations = c("***"),   
    vjust = 0.7,
    y_position = c(0.74),           
    tip_length = 0.05,
    textsize = 5
  ) +
  theme_minimal() +
  labs(title = "A) Sexuality", x = NULL, y = "Predicted Probability of Support", pattern = "Sexuality") +
  theme(
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5, margin = margin(10, 0, 0, 0)),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    axis.title.y = element_text(face="bold"),
    axis.text.x=element_blank()
    ) +
  ylim(y_limits)

plot2_data <- compute_summary(gt_national, "gt")

levels(plot2_data$gt) <- c("Conforming Man", "Slightly Nonconforming Man", "Nonconforming Man")

plot2 <- ggplot(plot2_data, aes(x = gt, y = mean_predicted_prob, fill = gt)) +
  geom_bar(stat = "identity", color = "black", width = 0.7) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, color = "black") + 
  geom_signif(
    comparisons = list(c("Conforming Man", "Nonconforming Man"), c("Slightly Nonconforming Man", "Nonconforming Man")),
    annotations = c("***", "***"),    
    vjust = 0.7,
    y_position = c(0.76, 0.70),           
    tip_length = 0.05,
    textsize = 5
  ) +
  scale_fill_grey(start = 0.3, end = 0.7) +  
  theme_minimal() +
  labs(title = "B) Gender Presentation", x = NULL, fill = "Gender Presentation") +
  theme(
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5, margin = margin(10, 0, 0, 0)),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x=element_blank(),
    axis.text.y=element_blank()
    ) +
  ylim(y_limits)

plot3_data <- compute_summary(gt_national, c("gt", "so"))

levels(plot3_data$so) <- c("Straight", "Gay")
levels(plot3_data$gt) <- c("Conforming Man", "Slightly Nonconforming Man", "Nonconforming Man")

plot3 <- ggplot(plot3_data, aes(x = gt, y = mean_predicted_prob, fill = gt, pattern = so)) +
  geom_bar_pattern(
    stat = "identity",
    position = position_dodge(width = 0.9), 
    color = "black",
    width = 0.7,
    pattern_fill = "black",
    pattern_density = 0.1,
    pattern_spacing = 0.05
    ) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, color = "black") +
  scale_fill_grey(start = 0.3, end = 0.7) +  
  theme_minimal() +
  labs(
    title = "C) Gender Presentation, By Sexuality",
    x = NULL,
    fill = "Gender Presentation",
    pattern = "Sexuality"
    ) +
  theme(
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5, margin = margin(10, 0, 0, 0)),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x=element_blank(),
    axis.text.y=element_blank(),
    strip.text = element_blank()
    ) +
  facet_wrap(~so, scales = "free_x") +  
  ylim(y_limits) +
  geom_signif(
    comparisons = list(c("Conforming Man", "Nonconforming Man")),
    annotations = c("***"),
    vjust = 0.7,
    y_position = c(0.76, 0.72),
    tip_length = 0.05,
    data = subset(plot3_data, so == "Straight")
  ) +
  geom_signif(
    comparisons = list(c("Conforming Man", "Nonconforming Man"), c("Slightly Nonconforming Man", "Nonconforming Man")),
    annotations = c("**", "***"),
    vjust = 0.7,
    y_position = c(0.73, 0.67),
    tip_length = 0.05,
    data = subset(plot3_data, so == "Gay")
  )

plot1_legend <- get_legend(plot1) 
plot2_legend <- get_legend(plot2) 

combined_legend <- plot_grid(plot1_legend, plot2_legend, ncol = 1, align = "v", rel_heights = c(1, 1))

plot1 <- plot1 + theme(legend.position = "none")
plot2 <- plot2 + theme(legend.position = "none")
plot3 <- plot3 + theme(legend.position = "none")

combined_plot <- plot_grid(
  plot_grid(plot1, plot2, plot3, ncol = 3, align = "h", rel_widths = c(1, 1, 2)),
  combined_legend,
  ncol = 2,
  rel_widths = c(3, 0.72)
)

combined_plot

#ggsave("Figure_4_March2025.pdf", combined_plot, width = 10, height = 5)

#ggsave("Figure_4_March2025.png", combined_plot, width = 10, height = 5)

```



```{r}

gt_national$predicted_prob <- predict(model_binary_national_1, type = "response")

mean_predicted_prob <- gt_national %>%
  group_by(gt) %>%
  summarise(mean_pred_prob = mean(predicted_prob, na.rm = TRUE))

# View the result
print(mean_predicted_prob)

#mass	0.6794921			(68%)
#neutral	0.6708913	(67%)
#fem	0.6138895	(61%)


mean_predicted_prob_so <- gt_national %>%
  group_by(so) %>%
  summarise(mean_predicted_prob_so = mean(predicted_prob, na.rm = TRUE))

print(mean_predicted_prob_so)

#straight	0.6874789			
#gay	0.6220131
```



#Predicted Probabilities PSSP

```{r}
# Remove rows with missing values in the DV or IVs used in the model
gt_pssp_clean <- gt_pssp %>%
  filter(!is.na(support_binary) & !is.na(so) & !is.na(gt) & !is.na(respondent) & !is.na(candidate))

# Fit the model with the cleaned data
model_binary_pssp_1 <- glmer(support_binary ~ so + gt + (1|respondent)+(1|candidate),
                             data = gt_pssp_clean, family = binomial(link = 'logit'))

# Predict probabilities
gt_pssp_clean$predicted_prob <- predict(model_binary_pssp_1, type = "response")

# Reorder levels of 'gt' (Gender Typicality) to be 'mass', 'neutral', 'fem'
gt_pssp_clean$gt <- factor(gt_pssp_clean$gt, levels = c("mass", "neutral", "fem"))

# Set common y-axis limits for all plots to ensure alignment (with a max of 0.8)
y_limits <- c(0, 0.8)

# 1. Plot 1: Predicted Probability by Sexuality using patterns
plot1_data_pssp_so <- gt_pssp_clean %>%
  group_by(so) %>%
  summarize(mean_predicted_prob = mean(predicted_prob))

plot1_data_pssp_gt <- gt_pssp_clean %>%
  group_by(gt) %>%
  summarize(mean_predicted_prob = mean(predicted_prob))

plot1_data_pssp_so
plot1_data_pssp_gt

#Sexuality (Straight: 0.6250671 vs. Gay: 0.7210530)
#Gender Typicality (0.6922143 vs. 0.6467443)

#0.04547 (0.6922143-0.6467443) 
```



```{r}
# New Dataframe (fem / ga as reference category)
gt_national_femref <- gt_national

gt_national_femref <- within(gt_national_femref, gt <- relevel(gt, ref = "fem"))
```

```{r}
#Checking significance for comparisons (Figure 4.C)

####################################################### All ##############################################################################

# All (Ref: GT)
model_binary_national_1 <- glmer(support_binary ~ so + gt + (1|respondent)+(1|candidate),
                              data = gt_national, family=binomial(link='logit'))

summary(model_binary_national_1)

#sogay       -0.44230    0.06380  -6.932 4.15e-12 ***
#gtneutral   -0.03888    0.07841  -0.496     0.62    
#gtfem       -0.30905    0.07807  -3.959 7.54e-05 ***

# All (Ref: GA)
model_binary_national_1_femref <- glmer(support_binary ~ so + gt + (1|respondent)+(1|candidate),
                              data = gt_national_femref, family=binomial(link='logit'))

summary(model_binary_national_1_femref)

#(Intercept)  1.08519    0.26431   4.106 4.03e-05 ***
#sogay       -0.44230    0.06380  -6.932 4.14e-12 ***
#gtmass       0.30905    0.07807   3.959 7.54e-05 ***
#gtneutral    0.27017    0.07778   3.473 0.000514 ***


#############################################################################  Gay ############################################################### 


# Gay (Ref: GT)
model_support_national_gt_gay <- glmer(support_binary ~ gt + (1|respondent),
                              data = subset(gt_national, so == "gay"), family=binomial(link='logit'))

summary(model_support_national_gt_gay)

#gtneutral    0.06372    0.13753   0.463  0.64316    
#gtfem       -0.42641    0.13674  -3.118  0.00182 ** 


# Gay (Ref: GA)

model_support_national_gt_gay_femref <- glmer(support_binary ~ gt + (1|respondent),
                              data = subset(gt_national_femref, so == "gay"), family=binomial(link='logit'))

summary(model_support_national_gt_gay_femref)

#gtmass        0.4264     0.1367   3.119 0.001817 ** 
#gtneutral     0.4901     0.1373   3.570 0.000357 ***

#############################################################################   Straight ################################################################# 

# Straight (Ref: GT)

model_support_national_gt_straight <- glmer(support_binary ~ gt + (1|respondent),
                              data = subset(gt_national, so == "straight"), family=binomial(link='logit'))

summary(model_support_national_gt_straight)

#gtneutral   -0.20176    0.11001  -1.834 0.066658 .  
#gtfem       -0.38863    0.11020  -3.527 0.000421 ***


# Straight (Ref: GA)

model_support_national_gt_straight_femref <- glmer(support_binary ~ gt + (1|respondent),
                              data = subset(gt_national_femref, so == "straight"), family=binomial(link='logit'))

summary(model_support_national_gt_straight_femref)

#gtmass       0.38864    0.11021   3.526 0.000421 ***
#gtneutral    0.18688    0.10856   1.721 0.085176 .  
```



# NOC & Anti-Gay Hostility as Moderators

```{r}
#Standardizing 0-1
range01 <- function(x,na.rm = T){(x-min(x, na.rm = T))/(max(x,na.rm = T)-min(x,na.rm = T))}

#NOC
gt_national$noc_standardized <- range01(gt_national$orcert.scale)
gt_pssp$noc_standardized <- range01(gt_pssp$noc.scale)

#Anti-Gay Attitudes
gt_national$homophobia_standardized <- range01(gt_national$homophobia.scale)
gt_pssp$homophobia_standardized <- range01(gt_pssp$homophobia.scale)


```


#Student Sample
```{r Student Sample, NOC and Homophobia}

#Support, NOC
supportModel_h3_gtsonoc_pssp <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + (1|respondent) + (1|candidate),
                              data = gt_pssp, family=binomial(link='logit')) 

summary(supportModel_h3_gtsonoc_pssp)

#Support, Homophobia 
supportModel_h4_sohom_pssp <- glmer(support_binary ~ so*homophobia_standardized + gt + (1|respondent) + (1|candidate),
                              data = gt_pssp, family=binomial(link='logit')) 

summary(supportModel_h4_sohom_pssp)

# Feeling, NOC

model_feeling_pssp_noc <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + (1 | respondent)+ (1 | candidate), 
                                           data = gt_pssp)
summary(model_feeling_pssp_noc)

# Feeling, Homophobia

model_feeling_pssp_hom <- lmer(feeling ~ so*homophobia_standardized + gt + (1 | respondent)+ (1 | candidate), 
                                           data = gt_pssp)
summary(model_feeling_pssp_hom)

```

#National Sample
```{r National Sample, NOC and Homophobia}

#Support, NOC

supportModel_h3_gtsonoc_national <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + (1|respondent) + (1|candidate),
                              data = gt_national, family=binomial(link='logit')) 

summary(supportModel_h3_gtsonoc_national)

#Support, Homophobia 

supportModel_h4_sohom_national <- glmer(support_binary ~ so*homophobia_standardized + gt + (1|respondent) + (1|candidate),
                              data = gt_national, family=binomial(link='logit')) 


summary(supportModel_h4_sohom_national)


# Feeling, NOC
model_feeling_national_noc <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + (1 | respondent)+ (1 | candidate), 
                                           data = gt_national)
summary(model_feeling_national_noc)

# Feeling, Homophobia

model_feeling_national_hom <- lmer(feeling ~ so*homophobia_standardized + gt + (1 | respondent)+ (1 | candidate), 
                                           data = gt_national)
summary(model_feeling_national_hom)


```

#Table 2: The Moderating Effect of NOC and Gay Attitudes on Biases Based on Sexuality and Gender Typicality

```{r}
# Apply the functions to the glmer support models
results_national_binary_noc <- extract_glmer_info(supportModel_h3_gtsonoc_national)
results_national_binary_hom <- extract_glmer_info(supportModel_h4_sohom_national)
results_pssp_binary_noc <- extract_glmer_info(supportModel_h3_gtsonoc_pssp)
results_pssp_binary_hom <- extract_glmer_info(supportModel_h4_sohom_pssp)


# Apply the function to the lmer feeling models
results_feeling_national_noc <- extract_lmer_info(model_feeling_national_noc)
results_feeling_national_hom <- extract_lmer_info(model_feeling_national_hom)
results_feeling_pssp_noc <- extract_lmer_info(model_feeling_pssp_noc)
results_feeling_pssp_hom <- extract_lmer_info(model_feeling_pssp_hom)


# Add model names to each dataframe
results_national_binary_noc$Model <- "National NOC (OR)"
results_national_binary_hom$Model <- "National Hom (OR)"
results_pssp_binary_noc$Model <- "PSSP NOC (OR)"
results_pssp_binary_hom$Model <- "PSSP Hom (OR)"
results_feeling_national_noc$Model <- "Feeling National NOC"
results_feeling_national_hom$Model <- "Feeling National Hom"
results_feeling_pssp_noc$Model <- "Feeling PSSP NOC"
results_feeling_pssp_hom$Model <- "Feeling PSSP Hom"


# Combine all results into one dataframe
combined_results_nochom <- rbind(results_pssp_binary_noc, results_pssp_binary_hom,results_feeling_pssp_noc, results_feeling_pssp_hom, results_national_binary_noc, results_national_binary_hom, results_feeling_national_noc, results_feeling_national_hom)

# Print combined results
print(combined_results_nochom)

```


```{r}

# Apply the format_model() custom function created earlier to the combined_results_nochom dataframe
combined_results_nochom <- combined_results_nochom %>%
  mutate(formatted_result = format_model(Estimate, Lower_95CI, Upper_95CI, P_value))

# Print the modified dataframe
print(combined_results_nochom)

# Reshape the dataframe to a wide format with the models as columns
wide_results_nochom <- combined_results_nochom %>%
  select(Term, Model, formatted_result) %>%
  tidyr::pivot_wider(names_from = Model, values_from = formatted_result)

# Print the reordered dataframe
print(wide_results_nochom)


```

```{r}
library(kableExtra)

# Create LaTeX table with the reordered columns
latex_table_nochom <- kable(wide_results_nochom, format = "latex", booktabs = TRUE, escape = FALSE,
                     col.names = c("Term", 
                                   "PSSP NOC (OR)", "PSSP Hom (OR)", 
                                   "Feeling PSSP NOC", "Feeling PSSP Hom", 
                                   "National NOC (OR)", "National Hom (OR)", 
                                   "Feeling National NOC", "Feeling National Hom")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down")) %>%
  add_header_above(c(" " = 1, "PSSP Models" = 4, "National Models" = 4))  # Custom header for the models


# Save as a LaTeX file
cat(latex_table_nochom, file = "table_2_nochom_results.tex")

```

#Second part of Table 2 (Variances and  group sizes)
```{r}
# Apply the extract_variances() custom function to all 8 NOC & Homophobia models
variance_pssp_binary_noc <- extract_variances_glmer(supportModel_h3_gtsonoc_pssp)
variance_pssp_binary_hom <- extract_variances_glmer(supportModel_h4_sohom_pssp)
variance_feeling_pssp_noc <- extract_variances_lmer(model_feeling_pssp_noc)
variance_feeling_pssp_hom <- extract_variances_lmer(model_feeling_pssp_hom)
variance_national_binary_noc <- extract_variances_glmer(supportModel_h3_gtsonoc_national)
variance_national_binary_hom <- extract_variances_glmer(supportModel_h4_sohom_national)
variance_feeling_national_noc <- extract_variances_lmer(model_feeling_national_noc)
variance_feeling_national_hom <- extract_variances_lmer(model_feeling_national_hom)


variance_pssp_binary_noc$Residual_Variance <- NA
variance_pssp_binary_hom$Residual_Variance <- NA
variance_national_binary_noc$Residual_Variance <- NA
variance_national_binary_hom$Residual_Variance <- NA


# Combine the results into a single dataframe
combined_variances_nochom <- rbind(variance_pssp_binary_noc,variance_pssp_binary_hom,variance_feeling_pssp_noc,variance_feeling_pssp_hom,variance_national_binary_noc,variance_national_binary_hom,variance_feeling_national_noc,variance_feeling_national_hom)

# Print combined variances dataframe
print(combined_variances_nochom)

```


```{r}
library(tidyr)
library(dplyr)

# Reshape the combined_variances_nochom dataframe so that the models are columns and variances are rows
reshaped_variances_nochom <- combined_variances_nochom %>%
  pivot_longer(cols = c(Respondent_Variance, Candidate_Variance, Residual_Variance, Num_Respondents, Num_Candidates), 
               names_to = "Metric", 
               values_to = "Value") %>%
  pivot_wider(names_from = Model, values_from = Value)

# Print reshaped variances
print(reshaped_variances_nochom)

```

```{r}
reshaped_variances_nochom <- reshaped_variances_nochom %>%
  mutate(across(where(is.numeric), ~ sprintf("%.2f", .)))


# Generate the LaTeX table
latex_table_variances_nochom <- reshaped_variances_nochom %>%
  kable(format = "latex", booktabs = TRUE, escape = FALSE, col.names = c("Metric", "PSSP Model 1 (OR)", "PSSP Model 2 (OR)", 
                                                                        "Feeling PSSP Model 1", "Feeling PSSP Model 2", 
                                                                        "National Model 1 (OR)", "National Model 2 (OR)", 
                                                                        "Feeling National Model 1", "Feeling National Model 2")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down"))

# Save the LaTeX table to a file (if you want)[[]]
cat(latex_table_variances_nochom, file = "reshaped_variances_table_nochom.tex")
```




#Partisanship as a Moderator

```{r}

#Student Sample

support_binary_pssp_party_3 <- glmer(support_binary ~ gt*partyid_3 + so*partyid_3 + (1|respondent) + (1|candidate), data = gt_pssp, family=binomial(link='logit'))
summary(support_binary_pssp_party_3)


feeling_pssp_party_3 <- lmer(feeling ~ gt*partyid_3 + so*partyid_3 + (1 | respondent)+ (1 | candidate), 
                                           data = gt_pssp)
summary(feeling_pssp_party_3)


#National Sample

support_binary_national_party_3 <- glmer(support_binary ~ gt*partyid_3 + so*partyid_3 + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))
summary(support_binary_national_party_3)


feeling_national_party_3 <- lmer(feeling ~ gt*partyid_3 + so*partyid_3 + (1 | respondent)+ (1 | candidate), data = gt_national)
summary(feeling_national_party_3)


```



#Table 3: The Effect of Candidate Sexuality and Gender Typicality on Electability, by Partisanship 

```{r}

# Apply the functions to the glmer support models
results_pssp_binary_partisanship <- extract_glmer_info(support_binary_pssp_party_3)
results_national_binary_partisanship <- extract_glmer_info(support_binary_national_party_3)

# Apply the function to the lmer feeling models
results_feeling_pssp_partisanship <- extract_lmer_info(feeling_pssp_party_3)
results_feeling_national_partisanship <- extract_lmer_info(feeling_national_party_3)


# Add model names to each dataframe
results_pssp_binary_partisanship$Model <- "PSSP Support"
results_feeling_pssp_partisanship$Model <- "PSSP Feeling"
results_national_binary_partisanship$Model <- "National Support"
results_feeling_national_partisanship$Model <- "National Feeling"




# Combine all results into one dataframe
combined_results_partisanship <- rbind(results_pssp_binary_partisanship,results_feeling_pssp_partisanship,results_national_binary_partisanship,results_feeling_national_partisanship)

# Print combined results
print(combined_results_partisanship)

```


```{r}

# Apply the format_model() custom function created earlier to the combined_results_nochom dataframe
combined_results_partisanship <- combined_results_partisanship %>%
  mutate(formatted_result = format_model(Estimate, Lower_95CI, Upper_95CI, P_value))

# Print the modified dataframe
print(combined_results_partisanship)

# Reshape the dataframe to a wide format with the models as columns
wide_results_partisanship <- combined_results_partisanship %>%
  select(Term, Model, formatted_result) %>%
  tidyr::pivot_wider(names_from = Model, values_from = formatted_result)

# Print the reordered dataframe
print(wide_results_partisanship)


```


```{r}
library(kableExtra)

# Create LaTeX table with the reordered columns
latex_table_partisanship <- kable(wide_results_partisanship, format = "latex", booktabs = TRUE, escape = FALSE,
                     col.names = c("Term", "PSSP Support", "PSSP Feeling", 
                                   "National Support", "National Feeling")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down")) %>%
  add_header_above(c(" " = 1, "PSSP Models" = 2, "National Models" = 2))  # Custom header for the models




# Save as a LaTeX file
cat(latex_table_partisanship, file = "table_3_partisanship_results.tex")

```

#Second part of Table 3 (Variances and  group sizes)
```{r}
# Apply the extract_variances() custom functions to all 4 partisanship moderators models
variance_pssp_support_partisanship <- extract_variances_glmer(support_binary_pssp_party_3)
variance_pssp_feeling_partisanship <- extract_variances_lmer(feeling_pssp_party_3)
variance_national_support_partisanship <- extract_variances_glmer(support_binary_national_party_3)
variance_national_feeling_partisanship <- extract_variances_lmer(feeling_national_party_3)


variance_pssp_support_partisanship$Residual_Variance <- NA
variance_national_support_partisanship$Residual_Variance <- NA


# Combine the results into a single dataframe
combined_variances_partisanship <- rbind(variance_pssp_support_partisanship,variance_pssp_feeling_partisanship,variance_national_support_partisanship,variance_national_feeling_partisanship)

# Print combined variances dataframe
print(combined_variances_partisanship)

```

```{r}
library(tidyr)
library(dplyr)

# Reshape the combined_variances_nochom dataframe so that the models are columns and variances are rows
reshaped_variances_partisanship <- combined_variances_partisanship %>%
  pivot_longer(cols = c(Respondent_Variance, Candidate_Variance, Residual_Variance, Num_Respondents, Num_Candidates), 
               names_to = "Metric", 
               values_to = "Value") %>%
  pivot_wider(names_from = Model, values_from = Value)

# Print reshaped variances
print(reshaped_variances_partisanship)

```

```{r}
reshaped_variances_partisanship <- reshaped_variances_partisanship %>%
  mutate(across(where(is.numeric), ~ sprintf("%.2f", .)))


# Generate the LaTeX table
latex_table_variances_partisanship <- reshaped_variances_partisanship %>%
  kable(format = "latex", booktabs = TRUE, escape = FALSE, col.names = c("Term", "PSSP Support", "PSSP Feeling", 
                                   "National Support", "National Feeling")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down"))

# Save the LaTeX table to a file (if you want)[[]]
cat(latex_table_variances_partisanship, file = "reshaped_variances_table_partisanship.tex")
```



# Subset Analyses by Partisanship (Analyses re: Figure 5)


```{r Feeling by Partisanship}

#National
feeling_national_Democrats <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), data = subset(gt_national, Democrat==1))
summary(feeling_national_Democrats)

feeling_national_Republicans <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), data = subset(gt_national, Republican==1))
summary(feeling_national_Republicans)

#PSSP

feeling_pssp_Democrats <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), data = subset(gt_pssp, partyid_3=="Democrat"))
summary(feeling_pssp_Democrats)

feeling_pssp_Republicans <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), data = subset(gt_pssp, partyid_3=="Republican"))
summary(feeling_pssp_Republicans)
```




```{r Support by Partisanship}

#National
support_binary_national_Democrats <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate),
                               data = subset(gt_national, Democrat==1), family=binomial(link='logit'))

summary(support_binary_national_Democrats)

support_binary_national_Republicans <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate),
                               data = subset(gt_national, Republican==1), family=binomial(link='logit'))

summary(support_binary_national_Republicans)

#PSSP

support_binary_pssp_Democrats <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate), data = subset(gt_pssp, partyid_3=="Democrat"), family=binomial(link='logit'))
summary(support_binary_pssp_Democrats)


support_binary_pssp_Republicans <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate), data = subset(gt_pssp, partyid_3=="Republican"), family=binomial(link='logit'))
summary(support_binary_pssp_Republicans)

```


```{r}
results_support_national_Democrats <- extract_glmer_info(support_binary_national_Democrats)
results_support_national_Republicans <- extract_glmer_info(support_binary_national_Republicans)
results_support_pssp_Democrats <- extract_glmer_info(support_binary_pssp_Democrats)
results_support_pssp_Republicans <- extract_glmer_info(support_binary_pssp_Republicans)


results_feeling_national_Democrats <- extract_lmer_info(feeling_national_Democrats)
results_feeling_national_Republicans <- extract_lmer_info(feeling_national_Republicans)
results_feeling_pssp_Democrats <- extract_lmer_info(feeling_pssp_Democrats)
results_feeling_pssp_Republicans <- extract_lmer_info(feeling_pssp_Republicans)


#Add Group Name (Student/National) to each dataframe
results_support_pssp_Democrats$group <- "Student Sample"
results_support_pssp_Republicans$group <- "Student Sample"
results_feeling_pssp_Democrats$group <- "Student Sample"
results_feeling_pssp_Republicans$group <- "Student Sample"

results_support_national_Democrats$group <- "National Sample"
results_support_national_Republicans$group <- "National Sample"
results_feeling_national_Democrats$group <- "National Sample"
results_feeling_national_Republicans$group <- "National Sample"

#Add Partisanship & DV to each dataframe

results_support_pssp_Democrats$partisanship <- "Democrats (DV: Support)"
results_support_pssp_Republicans$partisanship <- "Republicans (DV: Support)"
results_feeling_pssp_Democrats$partisanship <- "Democrats (DV: Feeling)"
results_feeling_pssp_Republicans$partisanship <- "Republicans (DV: Feeling)"

results_support_national_Democrats$partisanship <- "Democrats (DV: Support)"
results_support_national_Republicans$partisanship <- "Republicans (DV: Support)"
results_feeling_national_Democrats$partisanship <- "Democrats (DV: Feeling)"
results_feeling_national_Republicans$partisanship <-"Republicans (DV: Feeling)"


# Combine all results into one dataframe
combined_results_subset_partisanship <- rbind(results_support_national_Democrats,results_support_national_Republicans,results_support_pssp_Democrats,results_support_pssp_Republicans,results_feeling_national_Democrats,results_feeling_national_Republicans,results_feeling_pssp_Democrats,results_feeling_pssp_Republicans)

# Print combined results
print(combined_results_subset_partisanship)

plotdf_bothparties <- combined_results_subset_partisanship %>%
  filter(Term != "(Intercept)")

```




```{r}
library(ggplot2)
library(dplyr)
library(cowplot)
library(ggpubr) 

# Convert relevant variables to factors
plotdf_bothparties$partisanshipfactor <- as.factor(plotdf_bothparties$partisanship)
plotdf_bothparties$namesfactor <- as.factor(plotdf_bothparties$Term)

# Recode factor labels for clarity
plotdf_bothparties$namesfactor <- recode(plotdf_bothparties$namesfactor, 
                                         "sogay" = "Sexual Orientation \n (Gay)", 
                                         "gtneutral" = "Gender Presentation \n (Slightly Nonconforming Man)",
                                         "gtfem" = "Gender Presentation \n (Nonconforming Man)")

# Split data into two subsets
plotdf_support <- plotdf_bothparties %>% filter(grepl("Support", partisanshipfactor))
plotdf_feeling <- plotdf_bothparties %>% filter(grepl("Feeling", partisanshipfactor))

# Define base plot function
base_plot <- function(data, x_intercept) {
  ggplot(data, aes(x = Estimate, y = namesfactor, group = group, color = group, shape = group)) + 
    geom_point(size = 1.5, position = position_dodge(width = 0.40)) + 
    geom_errorbar(aes(xmin = Lower_95CI, xmax = Upper_95CI), size = 1, width = 0.2, position = position_dodge(width = 0.40)) + 
    geom_vline(xintercept = x_intercept, lty = 2) + 
    theme_bw() + 
    theme(plot.title = element_text(hjust = 0.4), 
          plot.subtitle = element_text(hjust = 0.4), 
          axis.text.y = element_text(color = "grey20", size = 11, face = "plain"), 
          legend.text = element_text(size = 11), 
          strip.text = element_text(size = 13),
          legend.position = "bottom",  # Ensure legend is placed at the bottom
          legend.box = "horizontal") +  # Force horizontal legend format
    xlab("") + 
    ylab("") + 
    facet_wrap(~partisanshipfactor, ncol = 2, scales = "free_x") + 
    labs(color = "Sample", shape = "Sample") + 
    scale_colour_grey(start = 0.30, end = 0.85)
}

# Create the two plots
plot_support <- base_plot(plotdf_support, x_intercept = 1) +
  scale_x_continuous(breaks = seq(0, 3, by = 0.5))  # Custom x-axis for support

plot_feeling <- base_plot(plotdf_feeling, x_intercept = 0)  # Default x-axis for feeling

# Combine plots using ggpubr::ggarrange()
combined_plot_bp <- ggarrange(plot_support, plot_feeling, 
                              ncol = 1, 
                              common.legend = TRUE,  # **Forces a single legend**
                              legend = "right")      # Puts legend at the bottom

# Display the final combined plot
print(combined_plot_bp)

# Save the plot
ggsave("plot_gtsoparties_gray_march2025.png", combined_plot_bp, width = 10, height = 7, dpi = 300, units = "in")

```



```{r}
# Create new dataframes based on conditions
gt_national_Democrats <- gt_national[gt_national$Democrat == 1, ]
gt_national_Republicans <- gt_national[gt_national$Republican == 1, ]

support_binary_national_Democrats <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate),
                               data = gt_national_Democrats, family=binomial(link='logit'))

summary(support_binary_national_Democrats)

support_binary_national_Republicans <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate),
                               data = gt_national_Republicans, family=binomial(link='logit'))

summary(support_binary_national_Republicans)

# Democrats

gt_national_Democrats$predicted_prob <- predict(support_binary_national_Democrats, type = "response")

mean_predicted_prob_so_Democrats <- gt_national_Democrats %>%
  group_by(so) %>%
  summarise(mean_predicted_prob_so_Democrats = mean(predicted_prob, na.rm = TRUE))

print(mean_predicted_prob_so_Democrats)

#straight	0.7222384	(72)
#gay	0.7793390	(78)
#Difference: 0.0571006 (6)

# Republicans

gt_national_Republicans$predicted_prob <- predict(support_binary_national_Republicans, type = "response")

mean_predicted_prob_so_Republicans <- gt_national_Republicans %>%
  group_by(so) %>%
  summarise(mean_predicted_prob_so_Republicans = mean(predicted_prob, na.rm = TRUE))

print(mean_predicted_prob_so_Republicans)

#straight	0.6992820			
#gay	0.4780086	
#Difference 0.2212734 (22)
```




##### SUPPLEMENTARY INFORMATION (SI) SECTION 6: REGRESSION RESULTS FOR CANDIDATE SUPPORT (1-6) ##### 


## Table SI-6.1 (H1 and H2 USING SUPPORT)

```{r}

supportModel_pssp <- ordinal::clmm(Support ~ so + gt + (1|respondent)+(1|candidate),
                              data = gt_pssp, link='logit')
summary(supportModel_pssp)

supportModel_h2_pssp <- ordinal::clmm(Support ~ so*gt + (1|respondent)+(1|candidate),
                              data = gt_pssp, link='logit')

summary(supportModel_h2_pssp)


supportModel_national <- ordinal::clmm(Support ~ so + gt + (1|respondent)+(1|candidate),
                              data = gt_national, link='logit')
summary(supportModel_national)

supportModel_h2_national <- ordinal::clmm(Support ~ so*gt + (1|respondent)+(1|candidate),
                              data = gt_national, link='logit')

summary(supportModel_h2_national)
```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(supportModel_pssp,supportModel_h2_pssp, supportModel_national, supportModel_h2_national),stars = c(0.001, 0.01, 0.05))

```


## Table SI-6.2 (Effect Moderation by NOC & Anti-Gay Hostility)

```{r}
#Standardizing 0-1
range01 <- function(x,na.rm = T){(x-min(x, na.rm = T))/(max(x,na.rm = T)-min(x,na.rm = T))}

#NOC
gt_national$noc_standardized <- range01(gt_national$orcert.scale)
gt_pssp$noc_standardized <- range01(gt_pssp$noc.scale)

#Anti-Gay Attitudes
gt_national$homophobia_standardized <- range01(gt_national$homophobia.scale)
gt_pssp$homophobia_standardized <- range01(gt_pssp$homophobia.scale)


```

```{r}

#NOC & Anti-Gay

supportModel_h3_gtsonoc_pssp <- ordinal::clmm(Support ~ so*noc_standardized + gt*noc_standardized + (1|respondent)+(1|candidate), data = gt_pssp, link='logit')

supportModel_h4_sohom_pssp <- ordinal::clmm(Support ~ so*homophobia_standardized + gt + (1|respondent)+(1|candidate), data = gt_pssp, link='logit')


supportModel_h3_gtsonoc_national <- ordinal::clmm(Support ~ so*noc_standardized + gt*noc_standardized + (1|respondent)+(1|candidate), data = gt_national, link='logit')


supportModel_h4_sohom_national <- ordinal::clmm(Support ~ so*homophobia_standardized + gt + (1|respondent)+(1|candidate), data = gt_national, link='logit')

```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(supportModel_h3_gtsonoc_pssp,supportModel_h4_sohom_pssp,supportModel_h3_gtsonoc_national,supportModel_h4_sohom_national),stars = c(0.001, 0.01, 0.05))

```


## Table SI-6.3 (Effect Moderation by Partisanship)


```{r}
supportModel_pssp_party_5 <- ordinal::clmm(Support ~ so*partyid_3 + gt*partyid_3 + (1|respondent)+(1|candidate), data = gt_pssp, link='logit')


supportModel_national_party_5 <- ordinal::clmm(Support ~ so*partyid_3 + gt*partyid_3 + (1|respondent)+(1|candidate), data = gt_national, link='logit')

```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(supportModel_pssp_party_5,supportModel_national_party_5),stars = c(0.001, 0.01, 0.05))

```


##### SUPPLEMENTARY INFORMATION (SI) SECTION 7: REGRESSION ANALYSES UNDERLYING FIGURES 4 AND 5 ##### 

# SECTION 7.1

```{r}
# New Dataframe (fem / gender nonconforming as reference category)
gt_national_femref <- gt_national

gt_national_femref <- within(gt_national_femref, gt <- relevel(gt, ref = "fem"))
```

#Note: In subsetted analyses (by gay or by straight), I include a random intercept only by respondent, rather than by both respondent and candidate. This decision was made after observing that models with both random intercepts fail to converge when subsetting the data by candidate sexuality (e.g., focusing only on gay candidates). This approach still adequately controls for within-respondent correlations while allowing for robust estimation of the fixed effects.

# Table SI-7.1 (Ref: GT)


```{r}

model_binary_national_all <- glmer(support_binary ~ so + gt + (1|respondent)+(1|candidate), data = gt_national, family=binomial(link='logit'))

model_binary_national_straight <- glmer(support_binary ~ gt + (1|respondent), data = subset(gt_national, so=="straight"), family=binomial(link='logit'))

model_binary_national_gay <- glmer(support_binary ~ gt + (1|respondent), data = subset(gt_national, so=="gay"), family=binomial(link='logit'))

```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(model_binary_national_all,model_binary_national_straight,model_binary_national_gay),stars = c(0.001, 0.01, 0.05))

```


# Table SI-7.2 (Ref: GA)
```{r}

model_binary_national_all_femref <- glmer(support_binary ~ so + gt + (1|respondent)+(1|candidate), data = gt_national_femref, family=binomial(link='logit'))

model_binary_national_straight_femref <- glmer(support_binary ~ gt + (1|respondent), data = subset(gt_national_femref, so=="straight"), family=binomial(link='logit'))

model_binary_national_gay_femref <- glmer(support_binary ~ gt + (1|respondent), data = subset(gt_national_femref, so=="gay"), family=binomial(link='logit'))


```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(model_binary_national_all_femref,model_binary_national_straight_femref,model_binary_national_gay_femref),stars = c(0.001, 0.01, 0.05))

```


# Table SI-7.3 (Democrats)

```{r Democrats}

#Student 


support_binary_pssp_Democrats <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate), data = subset(gt_pssp, partyid_3=="Democrat"), family=binomial(link='logit'))
summary(support_binary_pssp_Democrats)


feeling_pssp_Democrats <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), data = subset(gt_pssp, partyid_3=="Democrat"))
summary(feeling_pssp_Democrats)

#National 

support_binary_national_Democrats <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate),
                               data = subset(gt_national, Democrat==1), family=binomial(link='logit'))

summary(support_binary_national_Democrats)


feeling_national_Democrats <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), data = subset(gt_national, Democrat==1))
summary(feeling_national_Democrats)

```



```{r}
results_support_pssp_Democrats <- extract_glmer_info(support_binary_pssp_Democrats)
results_feeling_pssp_Democrats <- extract_lmer_info(feeling_pssp_Democrats)
results_support_national_Democrats <- extract_glmer_info(support_binary_national_Democrats)
results_feeling_national_Democrats <- extract_lmer_info(feeling_national_Democrats)


# Add model names to each dataframe
results_support_pssp_Democrats$Model <- "Support PSSP (OR)"
results_feeling_pssp_Democrats$Model <- "Feeling PSSP"
results_support_national_Democrats$Model <- "Support National (OR)"
results_feeling_national_Democrats$Model <- "Feeling National"


# Combine all results into one dataframe
combined_results_Democrats <- rbind(results_support_pssp_Democrats,results_feeling_pssp_Democrats,results_support_national_Democrats,results_feeling_national_Democrats)

# Print combined results
print(combined_results_Democrats)

```



```{r}

# Apply the format model function to the combined_results_Democrats dataframe
combined_results_Democrats <- combined_results_Democrats %>%
  mutate(formatted_result = format_model(Estimate, Lower_95CI, Upper_95CI, P_value))

print(combined_results_Democrats)

# Reshape the dataframe to a wide format with the models as columns
wide_results_Democrats <- combined_results_Democrats %>%
  select(Term, Model, formatted_result) %>%
  tidyr::pivot_wider(names_from = Model, values_from = formatted_result)

# Reorder columns
wide_results_Democrats <- wide_results_Democrats %>%
  select(Term, 
         "Support PSSP (OR)", ,"Feeling PSSP","Support National (OR)","Feeling National")

print(wide_results_Democrats)

```


```{r}
library(kableExtra)

# Create LaTeX table with the reordered columns
latex_table_Democrats <- kable(wide_results_Democrats, format = "latex", booktabs = TRUE, escape = FALSE,
                     col.names = c("Term", "Support PSSP (OR)","Feeling PSSP","Support National (OR)","Feeling National")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down")) %>%
  add_header_above(c(" " = 1, "PSSP Models" = 2, "National Models" = 2))  # Custom header for the models

cat(latex_table_Democrats, file = "table_Democrats_appendix.tex")

```


#Second part of Table (Variances and  group sizes)

```{r}

# Apply the function to all 4 models
variance_support_pssp_Democrats <- extract_variances_glmer(support_binary_pssp_Democrats)
variance_feeling_pssp_Democrats <- extract_variances_lmer(feeling_pssp_Democrats)
variance_support_national_Democrats<- extract_variances_glmer(support_binary_national_Democrats)
variance_feeling_national_Democrats <- extract_variances_lmer(feeling_national_Democrats)


variance_support_pssp_Democrats$Residual_Variance <- NA
variance_support_national_Democrats$Residual_Variance <- NA


# Combine the results into a single dataframe
combined_variances_Democrats <- rbind(variance_support_pssp_Democrats,variance_feeling_pssp_Democrats,variance_support_national_Democrats,variance_feeling_national_Democrats)

# Print combined_variances_Democrats dataframe
print(combined_variances_Democrats)

# Reshape the combined_variances_Democrats dataframe so that the models are columns and variances are rows
reshaped_variances_Democrats <- combined_variances_Democrats %>%
  pivot_longer(cols = c(Respondent_Variance, Candidate_Variance, Residual_Variance, Num_Respondents, Num_Candidates), 
               names_to = "Metric", 
               values_to = "Value") %>%
  pivot_wider(names_from = Model, values_from = Value)

# Print reshaped variances
print(reshaped_variances_Democrats)


```

```{r}
reshaped_variances_Democrats <- reshaped_variances_Democrats %>%
  mutate(across(where(is.numeric), ~ sprintf("%.2f", .)))


# Generate the LaTeX table
latex_table_variances_Democrats <- reshaped_variances_Democrats %>%
  kable(format = "latex", booktabs = TRUE, escape = FALSE, col.names = c("Metric", "Support PSSP (OR)","Feeling PSSP","Support National (OR)","Feeling National")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down"))



# Save the LaTeX table to a file (if you want)
cat(latex_table_variances_Democrats, file = "reshaped_variances_table_Democrats.tex")


```




# Table SI-7.4 (Republicans)

```{r Republicans}

#Student

support_binary_pssp_Republicans <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate), data = subset(gt_pssp, partyid_3=="Republican"), family=binomial(link='logit'))
summary(support_binary_pssp_Republicans)


feeling_pssp_Republicans <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), data = subset(gt_pssp, partyid_3=="Republican"))
summary(feeling_pssp_Republicans)

#National


support_binary_national_Republicans <- glmer(support_binary ~ so + gt + (1|respondent) + (1|candidate),
                               data = subset(gt_national, Republican==1), family=binomial(link='logit'))

summary(support_binary_national_Republicans)

feeling_national_Republicans <- lmer(feeling ~ so + gt + (1 | respondent)+ (1 | candidate), data = subset(gt_national, Republican==1))
summary(feeling_national_Republicans)


```


```{r}
results_support_pssp_Republicans <- extract_glmer_info(support_binary_pssp_Republicans)
results_feeling_pssp_Republicans <- extract_lmer_info(feeling_pssp_Republicans)
results_support_national_Republicans <- extract_glmer_info(support_binary_national_Republicans)
results_feeling_national_Republicans <- extract_lmer_info(feeling_national_Republicans)


# Add model names to each dataframe
results_support_pssp_Republicans$Model <- "Support PSSP (OR)"
results_feeling_pssp_Republicans$Model <- "Feeling PSSP"
results_support_national_Republicans$Model <- "Support National (OR)"
results_feeling_national_Republicans$Model <- "Feeling National"


# Combine all results into one dataframe
combined_results_Republicans <- rbind(results_support_pssp_Republicans,results_feeling_pssp_Republicans,results_support_national_Republicans,results_feeling_national_Republicans)

# Print combined results
print(combined_results_Republicans)


# Apply the format model function to the combined_results_Republicans dataframe
combined_results_Republicans <- combined_results_Republicans %>%
  mutate(formatted_result = format_model(Estimate, Lower_95CI, Upper_95CI, P_value))

# Print the modified dataframe
print(combined_results_Republicans)

# Reshape the dataframe to a wide format with the models as columns
wide_results_Republicans <- combined_results_Republicans %>%
  select(Term, Model, formatted_result) %>%
  tidyr::pivot_wider(names_from = Model, values_from = formatted_result)

# Reorder the columns in the desired order
wide_results_Republicans <- wide_results_Republicans %>%
  select(Term, 
         "Support PSSP (OR)", ,"Feeling PSSP","Support National (OR)","Feeling National")

# Print the reordered dataframe
print(wide_results_Republicans)

```


```{r}
library(kableExtra)

# Create LaTeX table with the reordered columns
latex_table_Republicans <- kable(wide_results_Republicans, format = "latex", booktabs = TRUE, escape = FALSE,
                     col.names = c("Term", "Support PSSP (OR)","Feeling PSSP","Support National (OR)","Feeling National")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down")) %>%
  add_header_above(c(" " = 1, "PSSP Models" = 2, "National Models" = 2))  # Custom header for the models

# Save as a LaTeX file
cat(latex_table_Republicans, file = "table_Republicans_appendix.tex")

```


#Second part of Table (Variances and  group sizes)

```{r}

# Apply the function to all 4 models
variance_support_pssp_Republicans <- extract_variances_glmer(support_binary_pssp_Republicans)
variance_feeling_pssp_Republicans <- extract_variances_lmer(feeling_pssp_Republicans)
variance_support_national_Republicans<- extract_variances_glmer(support_binary_national_Republicans)
variance_feeling_national_Republicans <- extract_variances_lmer(feeling_national_Republicans)


variance_support_pssp_Republicans$Residual_Variance <- NA
variance_support_national_Republicans$Residual_Variance <- NA


# Combine the results into a single dataframe
combined_variances_Republicans <- rbind(variance_support_pssp_Republicans,variance_feeling_pssp_Republicans,variance_support_national_Republicans,variance_feeling_national_Republicans)

# Print combined_variances_Republicans dataframe
print(combined_variances_Republicans)

# Reshape the combined_variances_Republicans dataframe so that the models are columns and variances are rows
reshaped_variances_Republicans <- combined_variances_Republicans %>%
  pivot_longer(cols = c(Respondent_Variance, Candidate_Variance, Residual_Variance, Num_Respondents, Num_Candidates), 
               names_to = "Metric", 
               values_to = "Value") %>%
  pivot_wider(names_from = Model, values_from = Value)

# Print reshaped variances
print(reshaped_variances_Republicans)


reshaped_variances_Republicans <- reshaped_variances_Republicans %>%
  mutate(across(where(is.numeric), ~ sprintf("%.2f", .)))


# Generate the LaTeX table
latex_table_variances_Republicans <- reshaped_variances_Republicans %>%
  kable(format = "latex", booktabs = TRUE, escape = FALSE, col.names = c("Metric", "Support PSSP (OR)","Feeling PSSP","Support National (OR)","Feeling National")) %>%
  kable_styling(latex_options = c("hold_position", "scale_down"))


# Save the LaTeX table to a file (if you want)
cat(latex_table_variances_Republicans, file = "reshaped_variances_table_Republicans.tex")


```


##### SUPPLEMENTARY INFORMATION (SI) SECTION 8: ANALYSES BY TREATMENT TYPE ##### 

#Gender Typicality: Face vs. Voice 

```{r}


feeling_voice_vs_face <- lmer(feeling ~ gt*voice + (1|respondent)+(1|candidate), data = gt_national)
summary(feeling_voice_vs_face)
 

supportModel_voice_vs_face  <- glmer(support_binary ~ gt*voice + (1|respondent) + (1 | candidate),
                                 data = gt_national, family = binomial(link = 'logit'))

summary(supportModel_voice_vs_face)


```

```{r}
texreg(l = list(supportModel_voice_vs_face, feeling_voice_vs_face),stars = c(0.001, 0.01, 0.05))

```



#Sexuality: Husband/Wife vs. Boyfriend/Girlfriend

```{r}
#In national sample, I only have one candidate with marriage==0, as such I use regular lm & glm models (as opposed to multilevel ones)

gt_national$marriage <- ifelse(gt_national$candidate == 1, 0, 1)

hw_vs_bg_main_feeling <- lm(feeling ~ so*marriage,  data = gt_national)
summary(hw_vs_bg_main_feeling)

hw_vs_bg_main_support <- glm(support_binary ~ so*marriage, data = gt_national, family=binomial(link='logit'))
summary(hw_vs_bg_main_support)

```


# Table SI-8.2: Moderation Analyses of Sexuality Penalty with Husband/Wife Cue 

```{r}
#Partisanship
hw_partisanship_support <- glm(support_binary ~ so*partyid_3, data = subset(gt_national, marriage==1), family=binomial(link='logit'))
summary(hw_partisanship_support)

hw_partisanship_feeling <- lm(feeling ~ so*partyid_3,  data = subset(gt_national, marriage==1))
summary(hw_partisanship_feeling)

#NOC

hw_noc_support <- glm(support_binary ~ so*noc_standardized, data = subset(gt_national, marriage==1), family=binomial(link='logit'))
summary(hw_noc_support)

hw_noc_feeling <- lm(feeling ~ so*noc_standardized,  data = subset(gt_national, marriage==1))
summary(hw_noc_feeling)

#Anti-Gay Attitudes

hw_homophobia_support <- glm(support_binary ~ so*homophobia_standardized, data = subset(gt_national, marriage==1), family=binomial(link='logit'))
summary(hw_homophobia_support)

hw_homophobia_feeling <- lm(feeling ~ so*homophobia_standardized,  data = subset(gt_national, marriage==1))
summary(hw_homophobia_feeling)


```

```{r}
texreg(l = list(hw_partisanship_support, hw_partisanship_feeling, hw_noc_support, hw_noc_feeling, hw_homophobia_support, hw_homophobia_feeling),stars = c(0.001, 0.01, 0.05))

```



# # Table SI-8.3: Moderation Analyses of Sexuality Penalty with Boyfriend / Girlfriend Cue

```{r}

#Partisanship

bg_partisanship_support <- glm(support_binary ~ so*partyid_3, data = subset(gt_national, marriage==0), family=binomial(link='logit'))

bg_partisanship_feeling <- lm(feeling ~ so*partyid_3,  data = subset(gt_national, marriage==0))

#NOC

bg_noc_support <- glm(support_binary ~ so*noc_standardized, data = subset(gt_national, marriage==0), family=binomial(link='logit'))

bg_noc_feeling <- lm(feeling ~ so*noc_standardized,  data = subset(gt_national, marriage==0))

#Anti-Gay Attitudes


bg_homophobia_support <- glm(support_binary ~ so*homophobia_standardized, data = subset(gt_national, marriage==0), family=binomial(link='logit'))

bg_homophobia_feeling <- lm(feeling ~ so*homophobia_standardized,  data = subset(gt_national, marriage==0))

```



```{r}
texreg(l = list(bg_partisanship_support,bg_partisanship_feeling,bg_noc_support,bg_noc_feeling,bg_homophobia_support,bg_homophobia_feeling),stars = c(0.001, 0.01, 0.05))

```




#Subsetted by Partisanship 

```{r Boyfriend vs Girfriend}

#Republicans

bg_Republicans_support <- glm(support_binary ~ so, data = subset(gt_national, marriage==0 & Republican==1), family=binomial(link='logit'))
summary(bg_Republicans_support)

bg_Republicans_feeling <- lm(feeling ~ so,  data = subset(gt_national, marriage==0 & Republican==1))
summary(bg_Republicans_feeling)


#Democrats

bg_Democrats_support <- glm(support_binary ~ so, data = subset(gt_national, marriage==0 & Democrat==1), family=binomial(link='logit'))
summary(bg_Democrats_support)

bg_Democrats_feeling <- lm(feeling ~ so,  data = subset(gt_national, marriage==0 & Democrat==1))
summary(bg_Democrats_feeling)


```


```{r Husband vs Wife}

#Republicans

hw_Republicans_support <- glm(support_binary ~ so, data = subset(gt_national, marriage==1 & Republican==1), family=binomial(link='logit'))
summary(hw_Republicans_support)

hw_Republicans_feeling <- lm(feeling ~ so,  data = subset(gt_national, marriage==1 & Republican==1))
summary(hw_Republicans_feeling)


#Democrats

hw_Democrats_support <- glm(support_binary ~ so, data = subset(gt_national, marriage==1 & Democrat==1), family=binomial(link='logit'))
summary(hw_Democrats_support)

hw_Democrats_feeling <- lm(feeling ~ so,  data = subset(gt_national, marriage==1 & Democrat==1))
summary(hw_Democrats_feeling)


```



##### SUPPLEMENTARY INFORMATION (SI) SECTION 10: ADDITIONAL ANALYSES ##### 

### ALTERNATIVE ANALYSES USING OLS AND LOGISTIC REGRESSION

#Table SI-10-1
```{r Student}
#Support
model_clm_support_pssp_h1 <- glm(support_binary ~ so + gt, data = gt_pssp, family=binomial(link='logit'))

model_clm_support_pssp_h2 <- glm(support_binary ~ so*gt, data = gt_pssp, family=binomial(link='logit'))

#Feeling

model_lm_feeling_pssp_h1 <- lm(feeling ~ so + gt, data = gt_pssp)

model_lm_feeling_pssp_h2 <- lm(feeling ~ so*gt, data = gt_pssp)

```

```{r National}
library(ordinal)

#Support
model_clm_support_national_h1 <- glm(support_binary ~ so + gt, data = gt_national, family=binomial(link='logit'))

model_clm_support_national_h2 <- glm(support_binary ~ so*gt, data = gt_national, family=binomial(link='logit'))


#Feeling

model_lm_feeling_national_h1 <- lm(feeling ~ so + gt, data = gt_national)

model_lm_feeling_national_h2 <- lm(feeling ~ so*gt, data = gt_national)
```


```{r}
#Custom function to exract random effects and report 
library("ordinal")
library("texreg")


texreg(l = list(model_clm_support_pssp_h1,model_clm_support_pssp_h2,model_lm_feeling_pssp_h1,model_lm_feeling_pssp_h2, model_clm_support_national_h1, model_clm_support_national_h2, model_lm_feeling_national_h1, model_lm_feeling_national_h2),stars = c(0.001, 0.01, 0.05))

```


#Table SI-10-2

#NOC & Anti-Gay Attitudes Moderators

```{r Student}
#Support
model_clm_support_pssp_h3 <- glm(support_binary ~ so*noc_standardized + gt*noc_standardized, data = gt_pssp, family=binomial(link='logit'))

model_clm_support_pssp_h4 <- glm(support_binary ~ so*homophobia_standardized + gt, data = gt_pssp, family=binomial(link='logit'))

#Feeling 

model_lm_feeling_pssp_h3 <- lm(feeling ~ so*noc_standardized + gt*noc_standardized, data = gt_pssp)

model_lm_feeling_pssp_h4 <- lm(feeling ~ so*homophobia_standardized + gt, data = gt_pssp)

```
 

```{r National}
library(ordinal)

#Support
model_clm_support_national_h3 <- glm(support_binary ~ so*noc_standardized + gt*noc_standardized, data = gt_national, family=binomial(link='logit'))

model_clm_support_national_h4 <- glm(support_binary ~ so*homophobia_standardized + gt, data = gt_national, family=binomial(link='logit'))


#Feeling 

model_lm_feeling_national_h3 <- lm(feeling ~ so*noc_standardized + gt*noc_standardized, data = gt_national)

model_lm_feeling_national_h4 <- lm(feeling ~ so*homophobia_standardized + gt, data = gt_national)
```


```{r}
library("ordinal")
library("texreg")

texreg(l = list(model_clm_support_pssp_h3, model_clm_support_pssp_h4, model_lm_feeling_pssp_h3, model_lm_feeling_pssp_h4, model_clm_support_national_h3, model_clm_support_national_h4, model_lm_feeling_national_h3, model_lm_feeling_national_h4),stars = c(0.001, 0.01, 0.05))

```

#Table SI-10-3

#Partisanship as Moderator

```{r Student}
model_clm_support_pssp_h5 <- glm(support_binary ~ so*partyid_3 + gt*partyid_3, data = gt_pssp, family=binomial(link='logit'))

model_lm_feeling_pssp_h5 <- lm(feeling ~ so*partyid_3 + gt*partyid_3, data = gt_pssp)

```

```{r National}
library(ordinal)

model_clm_support_national_h5 <- glm(support_binary ~ so*partyid_3 + gt*partyid_3, data = gt_national, family=binomial(link='logit'))

model_lm_feeling_national_h5 <- lm(feeling ~ so*partyid_3 + gt*partyid_3, data = gt_national)

```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(model_clm_support_pssp_h5,model_lm_feeling_pssp_h5,model_clm_support_national_h5,model_lm_feeling_national_h5),stars = c(0.001, 0.01, 0.05))


```


#Testing for 3-Way Interactions

#Table SI-10-4

```{r}
library(lme4)

#National

model_3_gtsopartisanship_feeling_national <- lmer(feeling ~ so*gt*partyid_3 + (1 | respondent)+ (1 | candidate), data = gt_national)
summary(model_3_gtsopartisanship_feeling_national)

model_3_gtsopartisanship_support_national <- glmer(support_binary ~ so*gt*partyid_3 + (1|respondent)+(1|candidate), data = gt_national, family=binomial(link='logit'))
summary(model_3_gtsopartisanship_support_national)

#Student
model_3_gtsopartisanship_feeling_student <- lmer(feeling ~ so*gt*partyid_3 + (1 | respondent)+ (1 | candidate), data = gt_pssp)
summary(model_3_gtsopartisanship_feeling_student)

model_3_gtsopartisanship_support_student <- glmer(support_binary ~ so*gt*partyid_3 + (1|respondent)+(1|candidate), data = gt_pssp, family=binomial(link='logit'))
summary(model_3_gtsopartisanship_support_student)

```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(model_3_gtsopartisanship_support_student,model_3_gtsopartisanship_feeling_student,model_3_gtsopartisanship_support_national,model_3_gtsopartisanship_feeling_national),stars = c(0.001, 0.01, 0.05))

```



##### SI SECTION 1: SAMPLE CHARACTERISTICS / SUMMARY STATS ##### 

#Table SI-1.1 (National Sample)
```{r}

#SEX
gt_national$sex_male <- NA
gt_national$sex_female <- NA
gt_national$sex_nborother <- NA

gt_national$sex_male <- ifelse(gt_national$sex =="Male", 1, 0)
gt_national$sex_female <- ifelse(gt_national$sex == "Female", 1, 0)
gt_national$sex_nborother <- ifelse(gt_national$sex == "Non-binary" | gt_national$sex == "Other (specify)", 1, 0)

#RACE
gt_national$race_white <- NA
gt_national$race_black <- NA
gt_national$race_latinx <- NA
gt_national$race_asian <- NA
gt_national$race_other <- NA

gt_national$race_white <- ifelse(gt_national$race_5 =="White", 1, 0)
gt_national$race_black <- ifelse(gt_national$race_5 == "Black", 1, 0)
gt_national$race_latinx <- ifelse(gt_national$race_5 == "Latinx", 1, 0)
gt_national$race_asian <- ifelse(gt_national$race_5 == "Asian", 1, 0)
gt_national$race_other <- ifelse(gt_national$race_5 == "Other", 1, 0)


#AGE (Binary)

gt_national$age_18_24 <- NA
gt_national$age_25_34 <- NA
gt_national$age_35_44 <- NA
gt_national$age_45_54 <- NA
gt_national$age_55_64 <- NA
gt_national$age_65_over <- NA

gt_national$age_18_24 <- ifelse(gt_national$age == 1, 1, 0)
gt_national$age_25_34 <- ifelse(gt_national$age == 2, 1, 0)
gt_national$age_35_44 <- ifelse(gt_national$age == 3, 1, 0)
gt_national$age_45_54 <- ifelse(gt_national$age == 4, 1, 0)
gt_national$age_55_64 <- ifelse(gt_national$age == 5, 1, 0)
gt_national$age_65_over <- ifelse(gt_national$age == 6, 1, 0)

#EDUCATION Binary

gt_national$edu_no_high_school <- NA
gt_national$edu_high_school <- NA
gt_national$edu_some_college <- NA
gt_national$edu_BA <- NA
gt_national$edu_graduate_degree <- NA

gt_national$edu_no_high_school <- ifelse(gt_national$edu2 =="Less than high school", 1, 0)
gt_national$edu_high_school <- ifelse(gt_national$edu2 == "High school only", 1, 0)
gt_national$edu_some_college <- ifelse(gt_national$edu2 == "Some college", 1, 0)
gt_national$edu_BA <- ifelse(gt_national$edu2 == "BA", 1, 0)
gt_national$edu_graduate_degree <- ifelse(gt_national$edu2 == "Graduate Degree", 1, 0)

#INCOME 
gt_national$income_ten<-recode(gt_national$income, "Less than $10,000"="Less than $10,000", 
                       "$10,000 - $19,999"="$10,000 - $19,999", "$20,000 - $29,999"="$20,000 - $29,999","$30,000 - $39,999"="$30,000 - $39,999", "$40,000 - $49,999"="$40,000 - $49,999","$50,000 - $59,999"="$50,000 - $59,999","$60,000 - $69,999"="$60,000 - $79,999","$70,000 - $79,999"="$60,000 - $79,999","$80,000 - $89,999"="$80,000 - $99,999","$90,000 - $99,999"="$80,000 - $99,999", "$100,000 - $149,999"="$100,000 - $149,999", "More than $150,000"="More than $150,000")

#INCOME NUMERIC
gt_national$income_numeric<-recode(gt_national$income, "Less than $10,000"=1, 
                       "$10,000 - $19,999"=2, "$20,000 - $29,999"=3,"$30,000 - $39,999"=4, "$40,000 - $49,999"=5,"$50,000 - $59,999"=6, "$60,000 - $69,999"=7,"$70,000 - $79,999"=7,"$80,000 - $89,999"=8,"$90,000 - $99,999"=8, "$100,000 - $149,999"=9, "More than $150,000"=10)

#INCOME (Binary)

gt_national$income_1 <- NA
gt_national$income_2 <- NA
gt_national$income_3 <- NA
gt_national$income_4 <- NA
gt_national$income_5 <- NA
gt_national$income_6 <- NA
gt_national$income_7 <- NA
gt_national$income_8 <- NA
gt_national$income_9 <- NA
gt_national$income_10 <- NA


gt_national$income_1 <- ifelse(gt_national$income_ten =="Less than $10,000", 1, 0)
gt_national$income_2 <- ifelse(gt_national$income_ten == "$10,000 - $19,999", 1, 0)
gt_national$income_3 <- ifelse(gt_national$income_ten == "$20,000 - $29,999", 1, 0)
gt_national$income_4 <- ifelse(gt_national$income_ten == "$30,000 - $39,999", 1, 0)
gt_national$income_5 <- ifelse(gt_national$income_ten == "$40,000 - $49,999", 1, 0)
gt_national$income_6 <- ifelse(gt_national$income_ten == "$50,000 - $59,999", 1, 0)
gt_national$income_7 <- ifelse(gt_national$income_ten == "$60,000 - $79,999", 1, 0)
gt_national$income_8 <- ifelse(gt_national$income_ten == "$80,000 - $99,999", 1, 0)
gt_national$income_9 <- ifelse(gt_national$income_ten == "$100,000 - $149,999", 1, 0)
gt_national$income_10 <- ifelse(gt_national$income_ten == "More than $150,000", 1, 0)

#PARTISANSHIP (Binary)
gt_national$Democrat <- NA
gt_national$Republican <- NA
gt_national$Independent <- NA

gt_national$Democrat <- ifelse(gt_national$partyid_3 =="Democrat", 1, 0)
gt_national$Republican <- ifelse(gt_national$partyid_3 =="Republican", 1, 0)
gt_national$Independent <- ifelse(gt_national$partyid_3 =="Independent", 1, 0)

table(gt_national$Democrat)
table(gt_national$Republican)
table(gt_national$Independent)

```


```{r}
library(dplyr)
gt_national_sample_stats <- gt_national %>% select(sex_male, sex_female, sex_nborother, race_white,race_black,race_latinx, race_asian, race_other, age_18_24, age_25_34, age_35_44, age_45_54, age_55_64,age_65_over,edu_no_high_school,edu_high_school,edu_some_college,edu_BA,edu_graduate_degree,income_1, income_2, income_3, income_4, income_5, income_6, income_7, income_8, income_9, income_10, Democrat, Republican, Independent, noc_standardized, homophobia_standardized)


```


```{r}
library(stargazer)
stargazer(gt_national_sample_stats, omit.summary.stat = c("p25", "p75"))
#because N is total N (respondents x each respondent seeing 4 candidates; the N in the appendix table is N/4 from the N that stargazer spits out)
#7884/4 = 1971
#7868/4 = 1967
```

#Table SI-1.2 Sample Statistics (Student Sample)

```{r}
unique(gt_pssp$sex)
#SEX
gt_pssp$sex_male <- NA
gt_pssp$sex_female <- NA
gt_pssp$sex_nborother <- NA

gt_pssp$sex_male <- ifelse(gt_pssp$sex =="Male", 1, 0)
gt_pssp$sex_female <- ifelse(gt_pssp$sex == "Female", 1, 0)
gt_pssp$sex_nborother <- ifelse(gt_pssp$sex == "Other (specify)", 1, 0)

#RACE
gt_pssp$race_white <- NA
gt_pssp$race_black <- NA
gt_pssp$race_latinx <- NA
gt_pssp$race_asian <- NA
gt_pssp$race_other <- NA

gt_pssp$race_white <- ifelse(gt_pssp$race_5 =="White", 1, 0)
gt_pssp$race_black <- ifelse(gt_pssp$race_5 == "Black", 1, 0)
gt_pssp$race_latinx <- ifelse(gt_pssp$race_5 == "Latinx", 1, 0)
gt_pssp$race_asian <- ifelse(gt_pssp$race_5 == "Asian", 1, 0)
gt_pssp$race_other <- ifelse(gt_pssp$race_5 == "Other", 1, 0)


#AGE (Binary)

gt_pssp$age_18_24 <- NA
gt_pssp$age_25_34 <- NA
gt_pssp$age_over_35 <- NA

gt_pssp$age_18_24 <- ifelse(gt_pssp$age_3 == 1, 1, 0)
gt_pssp$age_25_34 <- ifelse(gt_pssp$age_3 == 2, 1, 0)
gt_pssp$age_over_35 <- ifelse(gt_pssp$age_3 == 3, 1, 0)


#PARTISANSHIP (Binary)
gt_pssp$Democrat <- NA
gt_pssp$Republican <- NA
gt_pssp$Independent <- NA

gt_pssp$Democrat <- ifelse(gt_pssp$partyid_3 =="Democrat", 1, 0)
gt_pssp$Republican <- ifelse(gt_pssp$partyid_3 =="Republican", 1, 0)
gt_pssp$Independent <- ifelse(gt_pssp$partyid_3 =="Independent", 1, 0)


```



```{r}
library(dplyr)
gt_npssp_sample_stats <- gt_pssp %>% select(sex_male, sex_female, sex_nborother, race_white,race_black,race_latinx, race_asian, race_other, age_18_24, age_25_34, age_over_35, Democrat, Republican, Independent, noc_standardized, homophobia_standardized)


```


```{r}
library(stargazer)
stargazer(gt_npssp_sample_stats, omit.summary.stat = c("p25", "p75"))
#because N is total N (respondents x each respondent seeing 8 candidates; the N in the appendix table is N/8 from the N that stargazer spits out)
#4928/8 = 616
#4920/8 = 615
```


# SI SECTION 9 : MODERATORS OF BIAS EXPRESSION


```{r}
gt_national$edu5 <- gt_national$edu
gt_national$edu5<-recode(gt_national$edu5, "Less than high school degree"=1, 
                       "High school graduate (high school diploma including GED)"=2, "Some college but no degree"=3,"Associate degree (2-year)"=3,"Bachelor's degree (4-year)"=4,"Master's degree"=5,"Doctoral degree"=5,"Professional degree (JD, MD)"=5)

gt_national$sex_f <- as.factor(gt_national$sex)

gt_national$sex_f<-recode(gt_national$sex_f, "Female"="Female", 
                       "Male"="Male", "Non-binary"="Non-binary or Other","Other (specify)"="Non-binary or Other")

gt_national <- within(gt_national, sex_f <- relevel(sex_f, ref = "Female"))


#SEX

gt_national$sex_male_b <- NA

gt_national$sex_male_b <- ifelse(gt_national$sex_f == "Male", 1,
                        ifelse(gt_national$sex_f == "Female", 0, NA))


gt_national$age_s <- range01(gt_national$age)
gt_national$income_numeric_s <- range01(gt_national$income_numeric)
gt_national$edu5_s <- range01(gt_national$edu5)

```


```{r}

#SO & AGE
summary(lmer(feeling ~ so*age_s + (1|respondent)+(1|candidate), data = gt_national))
#sogay:age_s  -11.5793     1.3130 6519.9327  -8.819  < 2e-16 ***

#SO & GENDER

summary(lmer(feeling ~ so*sex_male_b + (1|respondent)+(1|candidate), data = gt_national))
#sogay:sex_male_b   -5.2923     0.8772 6380.4223  -6.033  1.7e-09 ***

#SO & EDUCATION

summary(lmer(feeling ~ so*edu5_s + (1|respondent)+(1|candidate), data = gt_national))
#sogay:edu5_s    6.534      1.685 6512.123   3.878 0.000106 ***
```


```{r}
#GT & AGE
summary(lmer(feeling ~ gt*age_s + (1|respondent)+(1|candidate), data = gt_national))

#GT & GENDER

summary(lmer(feeling ~ gt*sex_male_b + (1|respondent)+(1|candidate), data = gt_national))

#GT & EDUCATION

summary(lmer(feeling ~ gt*edu5_s + (1|respondent)+(1|candidate), data = gt_national))

```

#Table SI-9.1: NOC Moderation Feeling (Including Controls) 
```{r}
noc_basic_feeling <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + (1|respondent)+(1|candidate), data = gt_national)

noc_sex_feeling <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + sex_b + (1|respondent)+(1|candidate), data = gt_national)

noc_sex_age_feeling <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s +(1|respondent)+(1|candidate), data = gt_national)

noc_sex_age_edu_feeling <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s +edu5_s + (1|respondent)+(1|candidate), data = gt_national)
  
noc_sex_age_edu_income_feeling <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s +edu5_s + income_numeric_s + (1|respondent)+(1|candidate), data = gt_national)


noc_sex_age_edu_income_race_feeling <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s +edu5_s + income_numeric_s + race_5 + (1|respondent)+(1|candidate), data = gt_national)

noc_sex_age_edu_income_race_partisanship_feeling <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s +edu5_s + income_numeric_s + race_5 + partyid_3 + (1|respondent)+(1|candidate), data = gt_national)

noc_sex_age_edu_income_race_partisanship_homophobia_feeling <- lmer(feeling ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s + edu5_s + income_numeric_s + race_5 + partyid_3 + homophobia_standardized + (1|respondent)+(1|candidate), data = gt_national)
```

```{r}
library("ordinal")
library("texreg")


texreg(l = list(noc_basic_feeling,noc_sex_feeling,noc_sex_age_feeling,noc_sex_age_edu_feeling,noc_sex_age_edu_income_feeling,noc_sex_age_edu_income_race_feeling,noc_sex_age_edu_income_race_partisanship_feeling,noc_sex_age_edu_income_race_partisanship_homophobia_feeling),stars = c(0.001, 0.01, 0.05))

```


#Table SI-9.2: Anti-Gay Hostility Moderation Feeling (Including Controls)

```{r}
homophobia_basic_feeling <- lmer(feeling ~ so*homophobia_standardized + gt + (1|respondent)+(1|candidate), data = gt_national)

homophobia_sex_feeling <- lmer(feeling ~ so*homophobia_standardized + gt + sex_b + (1|respondent)+(1|candidate), data = gt_national)

homophobia_sex_age_feeling <- lmer(feeling ~ so*homophobia_standardized + gt + sex_b + age_s + (1|respondent)+(1|candidate), data = gt_national)

homophobia_sex_age_edu_feeling <- lmer(feeling ~ so*homophobia_standardized + gt + sex_b + age_s + edu5_s + (1|respondent)+(1|candidate), data = gt_national)
  
homophobia_sex_age_edu_income_feeling <- lmer(feeling ~ so*homophobia_standardized + gt + sex_b + age_s + edu5_s + income_numeric_s + (1|respondent)+(1|candidate), data = gt_national)


homophobia_sex_age_edu_income_race_feeling <- lmer(feeling ~ so*homophobia_standardized + gt + sex_b + age_s +edu5_s + income_numeric_s + race_5 + (1|respondent)+(1|candidate), data = gt_national)

homophobia_sex_age_edu_income_race_partisanship_feeling <- lmer(feeling ~ so*homophobia_standardized + gt + sex_b + age_s +edu5_s + income_numeric_s + race_5 + partyid_3 + (1|respondent)+(1|candidate), data = gt_national)

homophobia_sex_age_edu_income_race_partisanship_noc_feeling <- lmer(feeling ~ so*homophobia_standardized + gt + sex_b + age_s + edu5_s + income_numeric_s + race_5 + partyid_3 + noc_standardized + (1|respondent)+(1|candidate), data = gt_national)
```

```{r}
library("ordinal")
library("texreg")


texreg(l = list(homophobia_basic_feeling,homophobia_sex_feeling,homophobia_sex_age_feeling,homophobia_sex_age_edu_feeling,homophobia_sex_age_edu_income_feeling,homophobia_sex_age_edu_income_race_feeling,homophobia_sex_age_edu_income_race_partisanship_feeling,homophobia_sex_age_edu_income_race_partisanship_noc_feeling),stars = c(0.001, 0.01, 0.05))

```


#Table SI-9.3: NOC Moderation Support (Including Controls)
```{r}
noc_basic <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

noc_sex <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + sex_b + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

noc_sex_age <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

noc_sex_age_edu <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s + edu5_s + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

noc_sex_age_edu_income <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s + edu5_s + income_numeric_s + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

noc_sex_age_edu_income_race <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s + edu5_s + income_numeric_s + race_5 + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

noc_sex_age_edu_income_race_partisanship <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s + edu5_s + income_numeric_s + race_5 + partyid_3 + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

noc_sex_age_edu_income_race_partisanship_homophobia <- glmer(support_binary ~ so*noc_standardized + gt*noc_standardized + sex_b + age_s + edu5_s + income_numeric_s + partyid_3 + race_5 + homophobia_standardized + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))


```


```{r}
library("ordinal")
library("texreg")


texreg(l = list(noc_basic,noc_sex, noc_sex_age,noc_sex_age_edu,noc_sex_age_edu_income,noc_sex_age_edu_income_race,noc_sex_age_edu_income_race_partisanship,noc_sex_age_edu_income_race_partisanship_homophobia),stars = c(0.001, 0.01, 0.05))

```




#Table SI-9.4: Anti-Gay Hostility Moderation Support (Including Controls)

```{r}
homophobia_basic <- glmer(support_binary ~ so*homophobia_standardized + gt + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

homophobia_sex <- glmer(support_binary ~ so*homophobia_standardized + gt + sex_b + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

homophobia_sex_age <- glmer(support_binary ~ so*homophobia_standardized + gt + sex_b + age_s + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))


homophobia_sex_age_edu <- glmer(support_binary ~ so*homophobia_standardized + gt + sex_b + age_s + edu5_s + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

homophobia_sex_age_edu_income <- glmer(support_binary ~ so*homophobia_standardized + gt + sex_b + age_s + edu5_s + income_numeric_s + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

homophobia_sex_age_edu_income_race <- glmer(support_binary ~ so*homophobia_standardized + gt + sex_b + age_s + edu5_s + income_numeric_s + race_5 + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

homophobia_sex_age_edu_income_race_partisanship <- glmer(support_binary ~ so*homophobia_standardized + gt + sex_b + age_s + edu5_s + income_numeric_s + race_5 + partyid_3 + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))

homophobia_sex_age_edu_income_race_partisanship_noc <- glmer(support_binary ~ so*homophobia_standardized + gt + sex_b + age_s + edu5_s + income_numeric_s + partyid_3 + race_5 + noc_standardized + (1|respondent) + (1|candidate), data = gt_national, family=binomial(link='logit'))


```

```{r}
library("ordinal")
library("texreg")


texreg(l = list(homophobia_basic,homophobia_sex, homophobia_sex_age,homophobia_sex_age_edu,homophobia_sex_age_edu_income,homophobia_sex_age_edu_income_race,homophobia_sex_age_edu_income_race_partisanship,homophobia_sex_age_edu_income_race_partisanship_noc),stars = c(0.001, 0.01, 0.05))

```



##### SUPPLEMENTARY INFORMATION (SI) SECTION 3: MANIPULATION CHECKS AND TREATMENT CONTEXTUALIZATION ##### 

#Table SI-3.1 

```{r}
setwd("/Users/martinnaunov/Desktop/Desktop - Martin’s MacBook Pro (2)/Sexuality and Gender Typicality 2024")

pssp_manipulation_check <- read.csv("pssp_mc.csv")

```


```{r manipulation check}
library(lme4)
pssp_manipulation_check$candidate<-as.factor(pssp_manipulation_check$candidate)
pssp_manipulation_check$gt<-as.factor(pssp_manipulation_check$gt)
pssp_manipulation_check <- within(pssp_manipulation_check, gt <- relevel(gt, ref = "mass"))

options(contrasts = rep ("contr.treatment", 2))

pssp_check <- lmer(masculinity~gt+ (1|candidate),data=pssp_manipulation_check)

pssp_check_v1 <- lmer(masculinity~gt+ (1|candidate),data = subset(pssp_manipulation_check, v==1))

pssp_check_v0 <- lmer(masculinity~gt+ (1|candidate),data = subset(pssp_manipulation_check, v==0))

```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(pssp_check,pssp_check_v1,pssp_check_v0),stars = c(0.001, 0.01, 0.05))

```

#Means
```{r}
#Code to determine #'s in main text "The average rating shifted from 3.22 on a 7-point ....."

pssp_manipulation_factor <- pssp_manipulation_check %>%
  group_by(gt) %>%
  summarise(
    mean_gt = mean(masculinity, na.rm = TRUE),
    sd_gt = sd(masculinity, na.rm = TRUE)
  )

# Print the result
print(pssp_manipulation_factor)


```

#Table SI-3.2

# Cloud Research Pitch Manipulation Check (March 2024)

```{r}
setwd("/Users/martinnaunov/Desktop/Desktop - Martin’s MacBook Pro (2)/Sexuality and Gender Typicality 2024")

pitch_manipulation_check <- read.csv("cloud_research_pitch_manipulation_check.csv")

```


```{r}
pitch_manipulation_check$gt <- as.factor(pitch_manipulation_check$gt)

pitch_manipulation_check <- within(pitch_manipulation_check, gt <- relevel(gt, ref = "mass"))


mcheck_national_gt <- lmer(masculinity ~ gt + (1 | candidate), data = pitch_manipulation_check)
summary(mcheck_national_gt)

mcheck_national_emotional <- lmer(emotional_numeric ~ gt + (1 | candidate), data = pitch_manipulation_check)
summary(mcheck_national_emotional)

mcheck_national_ideo <- lmer(ideo ~ gt + (1 | candidate), data = pitch_manipulation_check)
summary(mcheck_national_ideo)

```

```{r}
library("ordinal")
library("texreg")

texreg(l = list(mcheck_national_gt,mcheck_national_emotional,mcheck_national_ideo),stars = c(0.001, 0.01, 0.05))

```


# Cloud Research (September 2024), Comparison with Real-World Gay Candidates

```{r}
setwd("/Users/martinnaunov/Desktop/Desktop - Martin’s MacBook Pro (2)/Sexuality and Gender Typicality 2024")

victory_fund <- read.csv("real_candidate_evals.csv")


```

```{r}
# Create a new dataframe by removing rows where real_life == 1 (i.e. dataframe containing only Victory Fund Headshots of real-world candidates)

real_candidates <- victory_fund %>% filter(real_candidate == 1)

# Create a new dataframe by removing rows where real_life == 0 (i.e. dataframe containing only my experimental stimuli)

stimuli <- victory_fund %>% filter(real_candidate == 0)

```


```{r}
library(dplyr)

#Victory Fund Headshots
average_masc_fem_victory <- real_candidates %>%
  group_by(CandidateID) %>%          # Group by CandidateID
  summarize(average_masc_fem = mean(masc_fem, na.rm = TRUE))  # Calculate the mean of masc_fem for each CandidateID

summary(average_masc_fem_victory$average_masc_fem)

#Mean: 3.096 (~3.1)

#Experimental Stimuli

average_masc_fem_stimuli <- stimuli %>%
  group_by(CandidateID) %>%          # Group by CandidateID
  summarize(average_masc_fem = mean(masc_fem, na.rm = TRUE))  # Calculate the mean of masc_fem for each CandidateID

summary(average_masc_fem_stimuli$average_masc_fem)
#3.084 (~3.1)

average_masc_fem_stimuli <- average_masc_fem_stimuli %>%
  mutate(gt = ifelse(CandidateID %% 2 == 0, "fem", "mass"))

```

```{r}
#Experimental Stimuli
manipulation_factor_stimuli <- average_masc_fem_stimuli %>%
  group_by(gt) %>%
  summarise(
    mean_gt = mean(average_masc_fem, na.rm = TRUE),
    sd_gt = sd(average_masc_fem, na.rm = TRUE)
  )

print(manipulation_factor_stimuli)

#fem	3.543900 (~3.5)
#mass	2.624188 (~2.6)

#Average femininity rating across the four "gender conforming" headshot stimuli used in my experiment was 2.6, while the "gender nonconforming" headshot stimuli had a mean rating of 3.5

```


#AVERAGES
```{r}
library(dplyr)

# Calculate the average masc_fem rating for each CandidateID and filter those with avg_masc_fem > 3.554467
candidates_above_355_threshold <- average_masc_fem_victory %>%
  filter(average_masc_fem > 3.543900)                               

# Count the number of candidates that meet the condition
num_candidates_above_355_threshold <- nrow(candidates_above_355_threshold)

(num_candidates_above_355_threshold/nrow(average_masc_fem_victory))*100
#24.5283 (25%)

#25% of gay candidates score above this on the masc/fem scale
#(Approximately 1 in 4 candidates score above this threshold)
#i.e. ~25% of gay candidates score more gender nonconforming than the average nonconforming stimulus version

```

```{r}
library(dplyr)

# Calculate the average masc_fem rating for each CandidateID and filter those with avg_masc_fem < 2.624188
candidates_below_262_threshold <- average_masc_fem_victory %>%
  filter(average_masc_fem < 2.624188)                               

num_candidates_below_262_threshold <- nrow(candidates_below_262_threshold)

(num_candidates_below_262_threshold/nrow(average_masc_fem_victory))*100
#31.44654% 

#About 31.4% of gay candidates are more masculine or gender conforming than the average gender conforming stimulus version)
#i.e. Almost 1 in 3 candidates score below this threshold

```



